use std::ops::Index;

pub fn min_edge_prime_num(number: u32) -> String {
    let mut k: u64 = 1;
    let mut prime_cnt: u64 = 0;
    let expected_res = number as f64 / 100.0;

    loop {
        k += 2;
        let delta = k - 1;
        let base = k * k;

        let mut i = 1;
        while i <= 3 {
            let candidate = base - i * delta;
            if miller_rabin(candidate) {
                prime_cnt += 1;
            }
            i += 1;
        }

        let res = prime_cnt as f64 / (2 * k - 1) as f64;
        if res < expected_res {
            return format!("{},{}", k, prime_cnt);
        }
    }
}

// Wrappers for witnesses
// Implementing Index trait without bound check
#[repr(transparent)]
#[derive(Debug, Clone, Copy, PartialEq, Eq, PartialOrd, Ord, Hash)]
pub struct Witnesses([(u64, Witness); 11]);

#[repr(transparent)]
#[derive(Debug, Clone, Copy, PartialEq, Eq, PartialOrd, Ord, Hash)]
pub struct Witness(&'static [u64]);

impl Index<usize> for Witnesses {
    type Output = (u64, Witness);

    fn index(&self, index: usize) -> &Self::Output {
        // no bound check
        unsafe {
            let ptr = self as *const Witnesses as *const (u64, Witness);
            &*ptr.add(index)
        }
    }
}

impl Index<usize> for Witness {
    type Output = u64;

    fn index(&self, index: usize) -> &Self::Output {
        // no bound check
        unsafe {
            let ptr = self.0.as_ptr();
            &*ptr.add(index)
        }
    }
}

// from crate: const-prime
// a high-performance miller-rabin prime test
pub fn miller_rabin(n: u64) -> bool {
    const NUM_BASES: usize = 11;
    /* const WITNESSES: [(u64, &[u64]); NUM_BASES] = [
        (2_046, &[2]),
        (1_373_652, &[2, 3]),
        (9_080_190, &[31, 73]),
        (25_326_000, &[2, 3, 5]),
        (4_759_123_140, &[2, 7, 61]),
        (1_112_004_669_632, &[2, 13, 23, 1_662_803]),
        (2_152_302_898_746, &[2, 3, 5, 7, 11]),
        (3_474_749_660_382, &[2, 3, 5, 7, 11, 13]),
        (341_550_071_728_320, &[2, 3, 5, 7, 11, 13, 17]),
        (3_825_123_056_546_413_050, &[2, 3, 5, 7, 11, 13, 17, 19, 23]),
        (u64::MAX, &[2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37]),
    ]; */
    const WITNESSES: Witnesses = Witnesses([
        (2_046, Witness(&[2])),
        (1_373_652, Witness(&[2, 3])),
        (9_080_190, Witness(&[31, 73])),
        (25_326_000, Witness(&[2, 3, 5])),
        (4_759_123_140, Witness(&[2, 7, 61])),
        (1_112_004_669_632, Witness(&[2, 13, 23, 1_662_803])),
        (2_152_302_898_746, Witness(&[2, 3, 5, 7, 11])),
        (3_474_749_660_382, Witness(&[2, 3, 5, 7, 11, 13])),
        (341_550_071_728_320, Witness(&[2, 3, 5, 7, 11, 13, 17])),
        (
            3_825_123_056_546_413_050,
            Witness(&[2, 3, 5, 7, 11, 13, 17, 19, 23]),
        ),
        (
            u64::MAX,
            Witness(&[2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37]),
        ),
    ]);

    if n == 2 || n == 3 {
        return true;
    } else if n <= 1 || n % 2 == 0 || n % 3 == 0 {
        return false;
    }

    // Use a small wheel to check up to log2(n).
    // This keeps the complexity at O(log(n)).
    let mut candidate_factor = 5;
    let trial_limit = n.ilog2() as u64;
    while candidate_factor <= trial_limit {
        if n % candidate_factor == 0 || n % (candidate_factor + 2) == 0 {
            return false;
        }
        candidate_factor += 6;
    }

    // Find r such that n = 2^d * r + 1 for some r >= 1
    let mut d = n - 1;
    while d % 2 == 0 {
        d >>= 1;
    }

    let mut i = 0;
    while i < NUM_BASES && WITNESSES[i].0 < n {
        i += 1;
    }
    let witnesses = WITNESSES[i].1;

    let mut i = 0;
    let len = witnesses.0.len();
    while i < len && witnesses[i] < n {
        if !miller_test(d, n, witnesses[i]) {
            return false;
        }
        i += 1;
    }

    true
}

/// Performs a Miller-Rabin test with the witness k.
const fn miller_test(mut d: u64, n: u64, k: u64) -> bool {
    let mut x = qpow(k, d, n);
    if x == 1 || x == n - 1 {
        return true;
    }

    while d != n - 1 {
        x = bmul(x, x, n);
        d *= 2;

        if x == 1 {
            return false;
        } else if x == n - 1 {
            return true;
        }
    }

    false
}

/// Calculates (a * b) mod m
const fn bmul(a: u64, b: u64, m: u64) -> u64 {
    (a as u128 * b as u128 % m as u128) as u64
}

/// Calculates (x^p) mod m
const fn qpow(mut x: u64, mut p: u64, m: u64) -> u64 {
    let mut ans = 1u64;
    x %= m;
    while p > 0 {
        if p % 2 == 1 {
            ans = bmul(ans, x, m);
        }
        x = bmul(x, x, m);
        p >>= 1;
    }
    ans
}
